model{
  for (i in start.week[1]:end.week[1]) {
    y[i] ~ dnorm(mu1[i], nu1[i])
    mu1[i] <- alpha1[weekid[i]] + gamma[survey.mode[i]] + delta[house[i]] + lambda[1] * cue.1[i] + lambda[2] * cue.2[i]
    nu1[i] <- pow((pow(tau[house[i]], design[i]) * sqrt(mu1[i] * (1 - mu1[i]) / n.obs[i])), -2)
  }
  for (i in start.week[2]:end.week[2]) {
    y[i] ~ dnorm(mu2[i], nu2[i])
    mu2[i] <- alpha2[weekid[i]] + gamma[survey.mode[i]] + delta[house[i]] + lambda[1] * cue.1[i] + lambda[2] * cue.2[i]
    nu2[i] <- pow((pow(tau[house[i]], design[i]) * sqrt(mu2[i] * (1 - mu2[i]) / n.obs[i])), -2)
  }
  for (i in start.week[3]:end.week[3]) {
    y[i] ~ dnorm(mu3[i], nu3[i])
    mu3[i] <- alpha3[weekid[i]] + gamma[survey.mode[i]] + delta[house[i]] + lambda[1] * cue.1[i] + lambda[2] * cue.2[i]
    nu3[i] <- pow((pow(tau[house[i]], design[i]) * sqrt(mu3[i] * (1 - mu3[i]) / n.obs[i])), -2)
  }
  for (i in start.week[4]:end.week[4]) {
    y[i] ~ dnorm(mu4[i], nu4[i])
    mu4[i] <- alpha4[weekid[i]] + gamma[survey.mode[i]] + delta[house[i]] + lambda[1] * cue.1[i] + lambda[2] * cue.2[i]
    nu4[i] <- pow((pow(tau[house[i]], design[i]) * sqrt(mu4[i] * (1 - mu4[i]) / n.obs[i])), -2)
  }
  for (i in start.week[5]:end.week[5]) {
    y[i] ~ dnorm(mu5[i], nu5[i])
    mu5[i] <- alpha5[weekid[i]] + gamma[survey.mode[i]] + delta[house[i]] + lambda[1] * cue.1[i] + lambda[2] * cue.2[i]
    nu5[i] <- pow((pow(tau[house[i]], design[i]) * sqrt(mu5[i] * (1 - mu5[i]) / n.obs[i])), -2)
  }
  for (i in start.week[6]:end.week[6]) {
    y[i] ~ dnorm(mu6[i], nu6[i])
    mu6[i] <- alpha6[weekid[i]] + gamma[survey.mode[i]] + delta[house[i]] + lambda[1] * cue.1[i] + lambda[2] * cue.2[i]
    nu6[i] <- pow((pow(tau[house[i]], design[i]) * sqrt(mu6[i] * (1 - mu6[i]) / n.obs[i])), -2)
  }
  for (i in start.week[7]:end.week[7]) {
    y[i] ~ dnorm(mu7[i], nu7[i])
    mu7[i] <- alpha7[weekid[i]] + gamma[survey.mode[i]] + delta[house[i]] + lambda[1] * cue.1[i] + lambda[2] * cue.2[i]
    nu7[i] <- pow((pow(tau[house[i]], design[i]) * sqrt(mu7[i] * (1 - mu7[i]) / n.obs[i])), -2)
  }
  for (i in start.week[8]:end.week[8]) {
    y[i] ~ dnorm(mu8[i], nu8[i])
    mu8[i] <- alpha8[weekid[i]] + gamma[survey.mode[i]] + delta[house[i]] + lambda[1] * cue.1[i] + lambda[2] * cue.2[i]
    nu8[i] <- pow((pow(tau[house[i]], design[i]) * sqrt(mu8[i] * (1 - mu8[i]) / n.obs[i])), -2)
  }
  alpha1[1] ~ dunif(lo[1], up[1])
  for (t in 2:n.weeks[1]) {
    alpha1[t] ~ dnorm(alpha1.star[t], omega.star[1])
    alpha1.star[t] <- alpha1[t - 1] + kappa * reshuffle[1, t]
  }
  alpha2[1] ~ dunif(lo[2], up[2])
  for (t in 2:n.weeks[2]) {
    alpha2[t] ~ dnorm(alpha2.star[t], omega.star[2])
    alpha2.star[t] <- alpha2[t - 1] + kappa * reshuffle[2, t]
  }
  alpha3[1] ~ dunif(lo[3], up[3])
  for (t in 2:n.weeks[3]) {
    alpha3[t] ~ dnorm(alpha3.star[t], omega.star[3])
    alpha3.star[t] <- alpha3[t - 1] + kappa * reshuffle[3, t]
  }
  alpha4[1] ~ dunif(lo[4], up[4])
  for (t in 2:n.weeks[4]) {
    alpha4[t] ~ dnorm(alpha4.star[t], omega.star[4])
    alpha4.star[t] <- alpha4[t - 1] + kappa * reshuffle[4, t]
  }
  alpha5[1] ~ dunif(lo[5], up[5])
  for (t in 2:n.weeks[5]) {
    alpha5[t] ~ dnorm(alpha5.star[t], omega.star[5])
    alpha5.star[t] <- alpha5[t - 1] + kappa * reshuffle[5, t]
  }
  alpha6[1] ~ dunif(lo[6], up[6])
  for (t in 2:n.weeks[6]) {
    alpha6[t] ~ dnorm(alpha6.star[t], omega.star[6])
    alpha6.star[t] <- alpha6[t - 1] + kappa * reshuffle[6, t]
  }
  alpha7[1] ~ dunif(lo[7], up[7])
  for (t in 2:n.weeks[7]) {
    alpha7[t] ~ dnorm(alpha7.star[t], omega.star[7])
    alpha7.star[t] <- alpha7[t - 1] + kappa * reshuffle[7, t]
  }
  alpha8[1] ~ dunif(lo[8], up[8])
  for (t in 2:n.weeks[8]) {
    alpha8[t] ~ dnorm(alpha8.star[t], omega.star[8])
    alpha8.star[t] <- alpha8[t - 1] + kappa * reshuffle[8, t]
  }
  gamma[1] <- gamma.star
  gamma[2] <- -gamma.star
  gamma.star ~ dnorm(0, 0.01)
  for (h in 1:10) {
    delta[h] <- delta.star[h]
  }
  delta[11] <- -delta.star[1] - delta.star[2] - delta.star[3] - delta.star[4] - 
    delta.star[5] - delta.star[6] - delta.star[7] - delta.star[8] - 
    delta.star[9] - delta.star[10]
  for (h in 1:10) {
    delta.star[h] ~ dnorm(0, 0.01)
  }
  for (p in 1:8) {
    omega.star[p] <- pow(omega[p], -2)
    omega[p] ~ dunif(0, 100)
  }
  for (h in 1:11) {
    tau[h] ~ dunif(0, 100)
  }
  kappa ~ dnorm(0, 0.01)
  lambda[1] ~ dnorm(0, 0.01)
  lambda[2] ~ dnorm(0, 0.01)
}